function SA = GIBBS1_SWR(S0,P0,P1,rho,T);

% function SA = GIBBS1_SWR(S0,P0,P1,rho,T);

% This file executes the Carter-Kohn backward sampler for the
% Stock-Watson-Romer model.

% S0,P0,P1 are outputs of the forward Kalman filter
A = [0 1 0; 0 1 0; 0 0 rho];

% initialize arrays for Gibbs sampler
SA = zeros(3,T); % artificial states
SM = zeros(3,1); % backward update for conditional mean of state vector
PM = zeros(3,3); % backward update for projection matrix
P = zeros(3,3); % backward update for conditional variance matrix
wa = randn(3,T); % draws for state innovations

% Backward recursions and sampling
% Terminal state
[V,D] = eig(P0(:,:,T));
Psqrt = real(V*(D.^.5)*V');
% if isreal(Psqrt) == 0,
%     'Warning: Matrix square root is complex'
% end
SA(:,T) = S0(:,T) + Psqrt*wa(:,T);
%SA(:,T) = S0(:,T) + real(sqrtm(P0(:,:,T)))*wa(:,T);

% iterating back through the rest of the sample
for i = 1:T-1,
   PM = (P0(:,:,T-i)*A')/P1(:,:,T-i);
   P = P0(:,:,T-i) - PM*A*P0(:,:,T-i);
   SM = S0(:,T-i) + PM*(SA(:,T-i+1) - A*S0(:,T-i));
   [V,D] = eig(P);
   Psqrt = real(V*(D.^.5)*V');
%    if isreal(Psqrt) == 0,
%     'Warning: Matrix square root is complex'
%    end
   SA(:,T-i) = SM + Psqrt*wa(:,T-i);
   %SA(:,T-i) = SM + real(sqrtm(P))*wa(:,T-i);
end